Variant Plateau’s law in atomically thin transition metal dichalcogenide dome networks

Since its fundamental inception from soap bubbles, Plateau’s law has sparked extensive research in equilibrated states. However, most studies primarily relied on liquids, foams or cellular structures, whereas its applicability has yet to be explored in nano-scale solid films. Here, we observed a variant Plateau’s law in networks of atomically thin domes made of solid two-dimensional (2D) transition metal dichalcogenides (TMDs). Discrete layer-dependent van der Waals (vdWs) interaction energies were experimentally and theoretically obtained for domes protruding in different TMD layers. Significant surface tension differences from layer-dependent vdWs interaction energies manifest in a variant of this fundamental law. The equivalent surface tension ranges from 2.4 to 3.6 N/m, around two orders of magnitude greater than conventional liquid films, enabling domes to sustain high gas pressure and exist in a fundamentally variant nature for several years. Our findings pave the way towards exploring variant discretised states with applications in opto-electro-mechanical devices.

shapes and an internal pressure of 80 atm. a, Perfectly round circular dome with a footprint radius = 500 nm results in a stiffness value (k) of 9.47 N/m b, An ellipse shape dome with 10% increment for the length in x-axis and 10% decrement for length in y-axis gives k = 9.27 N/m. c, An ellipse shape dome with 20% elongation of the length in x-axis and 20% decrement for length in y-axis gives k = 8.83 N/m. d, Dome with round rectangular edges, whose length in x axis is 1000 nm and length in y axis is 900nm gives k = 9.27 N/m.         4 . This is because in our work, the Poisson ratio of materials have been taken into account in the determination of Ah and Aτ, which affect the final results of E2D. Pressurised domes form due to the competition between inside pressure and interlayer vdWs energy, and this means dome samples are a good candidate to study the adhesion energy of layered materials, and the adhesion energy (γ) is given by 5,6 : where q is the dome profile factor where we use 2.1 obtained via statistical analysis, and ν is the materials' Poisson's ratio 5,6 . Using the stiffness values from layered dependent domes and Eq. (S1) to obtain the corresponding E2D, we calculated the γ in different layers of WS2 flakes to be 33.
where all variables have been previously introduced.
Using the equation above, we found that the pressure is exponential to 1/R and this trend is consistent with the fact that smaller domes have larger internal pressure. As shown in Supplementary Figure 4, different layered domes have distinguished capability to pressurise different amounts of hydrogen gas when their sizes are similar. 2L and 3L domes can compress far more hydrogen molecules within the same volume, leading to a difference in pressure of at least 50 atm compared with the same sized 1L dome. This trend is confirmed by numerical computation, and the results closely match the experimental values.

Strong boundary condition justification
In this work, a clamped boundary condition is applied as the 2D interface condition For the hemisphere (top) part, the volume can be described as: So the indentation volume (Vindentation) can be obtained by sum of two parts: Taking the dome sample shown in Fig S2 as  Comparing the indentation volume with dome's total volume, ( )/ , the indentation only creates less than 0.0002% volume deviation, which can be considered negligible. This is indicative of our reasoning for why we deliberately used smaller radius tips to reduce the influence of the indenters effect on domes during indentation. The van der Waals adhesion energy, thus, is believed to be quite capable of preventing slippage induced effects due to such small volume deviation. Given the AFM profiling over time in conjunction with our small indenter tip, this is highly supportive of a strongly clamped boundary condition.
Secondly, our boundary assumption is consistent with previous reports [5][6][7] . The clamped boundary condition is a universal boundary assumption in other works and our related previous work for both proton irradiation and stacking methods. In our previous work 5 , we used Raman to determine the Grunesian parameter and strain along the radial direction of domes that well matched the theoretical predictions, and the modelling assumption was based on a strongly clamped boundary condition.
Therefore, the clamped boundary condition employed in this work is valid, and the conclusions and calculations analytically and numerically support the experimental results with quite high accuracy.

Supplementary Information Note 2 Numerical Simulations using Finite Element Analysis (FEA)
For the COMSOL numerical implementation, COMSOL 5. approach. The identity edge pair location was determined from the experimental geometry profiles obtained via AFM. We simply had to use the measured geometry from AFM in terms of the domes boundaries as our COMSOL model uses finite element analysis, which only accounts for continuum mechanics and not a first principles energy calculation through molecular dynamics. This may be of interest for future work, however, the authors believe that the computational requirements would be extensive and angle minimisation would also need to be incorporated, therefore, COMSOL was a more practical approach for determining height profiles and stiffness of adjoining domes as shown throughout the manuscript and this supplemental work. From the comparisons in Fig. 4, and Supplementary Figures 8 and 9, it appears that this boundary condition was sufficient for numerically obtaining the height and stiffness profiles of coalescing domes. The same indentation technique was applied between the 2D and 3D models for consistency which in the low regime of stiffness measurements appears to match quite well with our experimental results.
We would also like to note, we specifically incorporated the use of a small indenter for the experiments which aided numerical simulations due to the indenter radius << dome footprint radius, thus heavily minimising nonlinear effects due to large surface area changes

Supplementary Information Note 3 Dome shape vs in-plane stiffness variation
To understand the influence of various dome footprints on measured stiffness, we performed a COMSOL simulations based on possible dome shapes as shown in Supplementary   Figure 3. The same amount of pressure is applied into the domes with varied geometries that deviate from a perfectly round dome such as ovals and a round-rectangular dome. The in-plane stiffness from the ideal dome with a circular footprint is 9.47 N/m, while the k values that have an oval shape with 10% and 20% shape deviation in the x-y directions have k values 9.27 and 8.83 N/m, respectively, indicating a 2.1% and 6.7% drop in stiffness. Moreover, the stiffness value from a round rectangular shape is also 9.27 N/m, leading to a 2.1% drop compared with the ideal dome scenario. Therefore, although the deviation of dome shape has some influence on the stiffness, the small deviation of in-plane stiffness indicates it is reasonable to use the round shape as a universal shape for the simulation and prediction of the in-plane stiffness from random domes. The circle in the centre is the indenter footprint and load location.

Supplementary Information Note 4 Pressure estimation via different approaches
To examine the presence of pressurised H2 molecules within the domes, the first method used is to place the H + irradiated WS2 dome sample in a low temperature environment (21 K) and monitor the inflation of domes with increased temperature. At low temperature (T = 21 K),  Comparing the pressure results governed by two different methods, they are comparable and have a good agreement, despite the small discrepancy which may arise from temperature induced modulated mechanical properties of TMD films. This indicates the inferred method employed in our work is capable to effectively probe the pressure condition in the domes.

Supplementary Table 2:
Capped gas pressure of selected domes estimated by two methods: hydrogen phase transition with temperature (blue) and nano-indentation method (green).

Supplementary Information Note 5 DFT simulations involving adhesion energy of different TMD layers
The plane-wave method in the framework of DFT as implemented in QUANTUM ESPRESSO code 10 is employed. Pseudopotentials are generated using the scheme proposed by Troullier and Martins (TM) 11 . GGA with the Perdew-Burke-Ernzerhof (PBE) 12

Supplementary Information Note 6 Bursting test to confirm the configuration of joint domes
To gather more direct evidence that a joint dome is made up of multiple domes within different basal planes of the TMD, a joint tri-dome with 1L, 2L and 3L dome regions

Mechanics of thin domes
In this section, a formulation for the mechanics of thin domes is presented. Domes are characterised by their height profile which has been shown to be universal for dome and tent like structures given by 5,6 ℎ = ℎ (1 − ( ) ) (S10) where h is the out of plane displacement, hm is the maximum height of the dome at the centre (as noted previously), r is the radial distance along a dome, R is the radius of the dome's footprint (also noted previously) and q is a constant. q has been shown to equal 2-2.2 for domes 5 (we use q = 2.1) and 2/3 for nanotents 6 . To solve the in-plane displacement field, the in plane equilibrium equation can be solved using the Foppl-von-Karman equation in the membrane limit given by 1,5 where u is the in-plane displacement. Whilst there have been several different formulations for the solution of Eq. (S11) depending on the assumptions of the out of plane displacement. 8,18 For this solution, we adopt the solution provided by Blundo et al. and Dai et al 5,6 . Solving out the second order differential equation for u, results in Having the radial displacement, the radial and circumferential strain fields along the radius of a dome can be formulated according to: Accordingly, the radial and circumferential strain fields in the strong shear region are formulated as 5,6 = ( , )ℎ 2 From the formulation of the strain field, the corresponding radial and circumferential stresses along the radius of a dome can be formulated as With the above definitions, a set of coupled nonlinear differential equations can be formulated for the solution of the transcendental equations given by Eq. (S11) and Eq. (S12)

Indentation of pressurised domes
In the following section, a discussion of the effects of indentation of the domes which are highly pressurised thin membranes is discussed. By realising that the effect of a finite radius indenter compared to the radius of the dome results in a logarithmic term 1  In this work, we fabricated 1L, 2L and 3L domes and it is expected that the 2D modulus should increase with layer number and this has been shown for suspended materials subjected to AFM indentation 20 . It is expected that E2D is larger for multiple layer domes resulting in increased pressure to achieve the same height and an increased stiffness (see Eqs. (S17), (S18), (S24) and (S25)). Moreover, for larger E2D in 2L and 3L systems, it is expected that that corresponding γ also increases due to the exfoliation energy required to separate multiple layers from the initial few basal plane layers.

Supplemental Information Note 8 Plateau's law in liquid bubbles and variations in joint angle in analogous nano dome networks
Energy minimisation problems exist everywhere in the universe from large scale systems, macro, micro and in nanoscale systems. Such energy minimisation topologies are abundant in nature in the formation of cell structures and bees forming honeycombs. Probably, the most famous energy minimisation problem is the formation of liquid soap bubbles and the interaction of merging bubbles to assume a minimum energy formation. As bubbles merge, smaller bubbles with larger pressure always protrude into larger bubbles. Similar to the hexagonal structure of graphene, liquid soap bubbles form at 120° angles to assume a minimum energy shape and this is considered the double bubble conjecture 21 . However, if soapy liquid merges with four vertices then they form the Maraldi angle at 109.5°. It is very easy to observe that liquid bubbles always rearrange themselves for a minimum energy state in terms of their surface area to volume ratio. From Laplace's law, we have the pressure inside a half sphere is given by: where Δp is the pressure differential between inner and outer wall of a bubble (same terminology for internal pressure for a dome), σ is the surface tension at the liquid surface interface (N/m) and Rs is the Radius of the sphere. where all previous terms are the same as previously introduced in Note 1 (Eqs. (S2-S5)) and After rigorous searches and statistics of nano domes in the literature, q was found to be within 2.0-2.2 5 , which is consistent with the our current work and a value of 2.1±0.05 was extracted using Eq. (S10) Eq. (S31) can be rearranged in a similar format analogous to Laplace's law for liquid bubbles given by Δ = 2 (S33) the σ is given by However, since our domes have negligible bending stiffness, the effective surface tension is given by For the surface tension of the bottom single layer on the shared boundary, the exfoliation process should mean that 1 ′ represents a portion of 1 as it should require less energy to detach a 1L sheet from a 2L sheet rather than from the bulk. As the shared boundary can be regarded as the exfoliation of a top single layer from a bilayer film, we propose the adhesion energy for this process is calculated as: where 1 and 2 are the adhesion energy of 1L and 2L domes, respectively. Similarly, the surface tension of the shared boundary can be calculated as It should also be noted that for the 3L dome in a situation of 3 merging bubbles appears to be returning to the normal scenario, therefore, future work with rigorous searching conditions and if domes in deeper layers can be generated should return to Plateau's law with joint angles of 120°.